import numpy as np
import matplotlib.pyplot as plt

# 参数设置
fs = 1  # 采样率 = 1 Hz (每秒1个点)
N = 100  # 总点数 = 100
pulse_interval = 10  # 脉冲间隔 = 10个点

# 1. 生成脉冲信号
t = np.arange(N)  # 时间序列 [0, 1, 2, ..., 99] 秒
signal = np.zeros(N)
signal[::pulse_interval] = 1  # 每10个点设置一个脉冲（索引0,10,20,...,90）

# 2. 计算FFT
fft_result = np.fft.fft(signal)  # 复数频谱
freq = np.fft.fftfreq(N, d=1/fs)  # 频率轴 (Hz)
magnitude = np.abs(fft_result) / N * 2  # 幅度谱归一化 <cite data-id='30008'>30008</cite>
magnitude[0] /= 2  # 直流分量特殊处理 

# 3. 二次FFT
fft_result2 = np.fft.fft(fft_result)  # 复数频谱
magnitude2 = np.abs(fft_result2) / N * 2  # 幅度谱归一化 <cite data-id='30008'>30008</cite>
magnitude2[0] /= 2  # 直流分量特殊处理 

# 3. 可视化
fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(12, 8))

# 时域信号
ax1.stem(t, signal, linefmt='b-', markerfmt='bo', basefmt=" ")
ax1.set_title("Periodic Pulse Signal (Time Domain)")
ax1.set_xlabel("Time (s)")
ax1.set_ylabel("Amplitude")
ax1.grid(alpha=0.3)
ax1.set_xticks(np.arange(0, N, 10))

# 频域幅度谱
positive_freq_mask = (freq >= 0) & (freq <= fs/2)  # 取正频率部分 (0~0.5 Hz)
ax2.plot(freq[positive_freq_mask], magnitude[positive_freq_mask], 'r-o')
ax2.set_title("FFT Magnitude Spectrum")
ax2.set_xlabel("Frequency (Hz)")
ax2.set_ylabel("Normalized Amplitude")
ax2.grid(alpha=0.3)
ax2.set_xticks(np.arange(0, 0.51, 0.1))

ax3.plot(freq[positive_freq_mask], magnitude2[positive_freq_mask], 'r-x')
ax3.set_title("FFTx2 Magnitude Spectrum")
ax3.set_xlabel("Frequency (Hz)")
ax3.set_ylabel("Normalized Amplitude - FFT*2")
ax3.grid(alpha=0.3)
ax3.set_xticks(np.arange(0, 0.51, 0.1))

plt.tight_layout()
plt.show()